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Practical applications of thermoacoustic tomography require numerical inversion of 

the spherical mean Radon transform with the centers of integration spheres occupying 

an open surface. Solution of this problem is needed (both in 2-D and 3-D) because 

"^ frequently the region of interest cannot be completely surrounded by the detectors, as 

^ it happens, for example, in breast imaging. We present an efficient numerical algorithm 

Ch for solving this problem in 2-D (similar methods are applicable in the 3-D case). Our 

method is based on the numerical approximation of plane waves by certain single layer 

potentials related to the acquisition geometry. After the densities of these potentials 

have been precomputed, each subsequent image reconstruction has the complexity of 

I/-] the regular filtration backprojection algorithm for the classical Radon transform. The 

C^ peformance of the method is demonstrated in several numerical examples: one can see 

that the algorithm produces very accurate reconstructions if the data are accurate and 

L^ sufficiently well sampled, on the other hand, it is sufficiently stable with respect to noise 

(yQ in the data. 

o 

> 
|/J Introduction 

C^ The thermoacoustic tomography (TAT) is based on measurements of acoustic waves excited 

in the patient's body by an external electromagnetic pulse [15,36]. The method is gaining 
popularity among researchers because it combines the high resolution of the ultrasound to- 
mography with the high contrast of the images attainable due to the strong variance in the 
electromagnetic properties of the body tissues. One of the most promising applications of 
this technique is breast imaging, where absorption of the electromagnetic energy in tumors 
is several times higher than that in healthy tissues. Under certain simplifying assumptions 
(most importantly, that of a constant speed of sound in the tissue) the corresponding recon- 
struction problem can be reduced to the inversion of the spherical mean Radon transform. In 
the present paper we propose an efficient numerical algorithm for solving this problem in 2-D 
in the case when the detectors lie on an open curve only partially surrounding the region of 
interest. 



A significant progress has been achieved recently in mathematics of TAT; an extended 
discussion and relevant references can be found in reviews [1,12,16]. In particular, several 
versions of explicit inversions formulas [8, 11, 12, 18, 37] and certain series solutions [19, 22, 23] 
have been obtained for measurement schemes with the detectors lying either on closed surfaces 
surrounding the region of interest (ROI), or on certain open but unbounded surfaces, such as 
a plane or infinite cylinder. 

However, in most applications the ROI is not the whole human body but rather a part of 
it — as it happens, for example, in mammography. In such situations only some part of the 
region boundary can be covered by the detectors. This, in turn, necessitates the development 
of algorithms that can reconstruct the image from such data. 

The simplest approach to image reconstruction from open surfaces is to place detectors on 
a truncated plane or cylinder and to apply the known inversion formula(s) for the full plane 
or infinite cylinder. Of course, the reconstruction would not be exact; moreover, it has been 
shown that for a stable reconstruction of a compactly supported function from its spherical 
means it is necessary to know the integrals over sufficiently many spheres, so that the so-called 
"visibility" condition is satisfied [20,24,25,38,39]. (It have been shown, in particular, that 
under this condition the reconstruction problem can be reduced either to the inversion of an 
elliptic operator or, after a certain filtration, to a solution of the Fredholm integral equation 
of the second kind [17, 24, 25, 27]. In both cases the problem is only mildly unstable, similarly 
to the inversion of the standard Radon transform [21]). The "visibility" condition (for TAT) 
is satisfied if for each point x in the ROI each straight line passing through x intersects the 
measuring surface at least once. A sphere surrounding the object satisfies this condition while 
an infinite plane and an infinite cylinder do not. For a truncated plane or truncated cylinder 
the set of "bad" directions is even larger than for their infinite counterparts. If this condition 
is violated, it is practically impossible to accurately reconstruct sharp interfaces (material 
boundaries) at those points of the image for which some of the straight lines normal to the 
interface do not intersect the measuring surface. (In 2-D an exact reconstruction technique 
was developed in [26] for the case when the centers of the integration circles lie on a segment 
of a straight line. However, the reconstruction is still, in general, unstable in this case.) 

For a given bounded ROI there exist many bounded open acquisition surfaces that do 
satisfy the visibility condition. Almost all known reconstruction techniques applicable to such 
surfaces are of approximate nature. For example, by "approximating" the integration spheres 
by planes and by applying some version of the classical inverse Radon transform, one can 
reconstruct an "approximation" to the image. Due to the symmetry in the classical Radon 
projections, the normals to the integration planes should fill only a half of a unit sphere, making 
possible the reconstruction from an open measurement surface. A more sophisticated approach 
is represented by the so-called "straightening" methods [33,34] based on the approximate 
reconstruction of the classical Radon projections from the measurements that correspond to 
the spherical mean Radon transform. However, these methods do not yield the exact solution 
but rather a parametrix of the problem; the image is reconstructed only up to a certain 
smooth term. In other words, while jumps corresponding to sharp material interfaces are 
reconstructed accurately, the accuracy of the lower spatial frequencies can not be guaranteed. 
Unlike the approximations resulting from the discretization of the exact inversion formulas 
(in the situations when such formulas are known), the parametrix approximations do not 
converge when the discretization of the data is refined (in the absence of noise). In [28], an 
exact inversion formula for the spherical surface is used to obtain approximate reconstructions 



from other measurement surfaces; in order to further improve the approximation, additional 
corrections are introduced. These methods yield different parametrices than the one proposed 
in [33,34]; which one is better remains an open question. 

An accurate numerical reconstruction from an open surface acquisition scheme was demon- 
strated in [4] (see also [31]). In this work a version of an iterative algebraic reconstruction 
algorithm was successfully employed to recover a numerically generated phantom. Iterative 
algebraic reconstruction algorithms are, however, notoriously slow; the above-mentioned re- 
construction, for example, required the use of a cluster of computers and took 100 iterations 
to converge. A faster converging algorithm can be obtained by combining iterative refining 
with a parametrix-type algorithm [28,30]. These methods are closely related to the general 
scheme proposed in [5] for the inversion of the generalized Radon transform with integration 
over general manifolds. It reduces the problem to the Fredholm integral equation of the second 
kind, which is well suited for numerical solution. Such an approach can be viewed as using a 
parametrix method as an efficient preconditioner for an iterative solver, which, as it has been 
shown in numerical experiments, significantly accelerates the convergence of the iterations.. 

Finally, in [32] an interesting attempt has been made to generate the absent data from the 
consistency conditions on the spherical mean Radon transform, in order to "numerically close" 
the open measurement surface. The resulting algorithm, however, seems to be less accurate 
than methods we mentioned earlier, and it exhibits instability on higher spatial frequencies. 

In the present paper we propose a novel non-iterative algorithm for the numerical inversion 
of the spherical mean Radon transform from open surfaces in 2-D. Our numerical experiments 
show that this method yields very accurate reconstructions in the absence of noise, and is 
almost as stable as the classical filtration/backprojection (FBP) algorithm (widely used for 
the inversion of the regular Radon transform [14,21]). The present algorithm requires pre- 
computation of certain functions that jointly serve as a numerical filter; this needs to be done 
only once for a given configuration of detectors. After these functions have been computed, 
each reconstruction has numerical complexity similar to that of the FBP. 

1 Formulation of the problem 

In the thermoacoustic tomography acoustic detectors measure the pressure of the outgoing 
wave radiating from the patient's body. This acoustic wave is generated by the thermoacoustic 
expansion of the tissues initiated by a very short electromagnetic pulse. The initial pressure 
f{x) strongly depends on the type of tissue, and is significantly higher for tumors, since they 
happen to absorb much more electromagnetic energy than healthy tissue. Thus, recovering 
f{x) would yield important medical information about the location and shape of tumors. 

Let us denote by g{z,t) the pressure registered at the moment t by the acoustic detector 
placed at the point z. Under certain assumptions (such as a constant sound speed in the body, 
ideal infinitely small detector, infinitely short pulse, and so on), one can recover from the 
measurements the integral of f{x) over a sphere of radius r = ct centered at z. A dimensional 
analysis shows that in order to reconstruct a function of a 3-D variable, the detectors should 
cover some measurement surface S. Assuming for simplicity that the speed of sound equals 
unity, one arrives at the following formulation of the inverse problem: reconstruct function 
f{x) supported within some 3-D region of interest Q from known values of its integrals g{z, r) 
over all spheres of radius r with centers z lying on surface S {S and Q are disjoint sets). 




Figure 1: Geometry 

A 2-D version of this problem also arises in the thermoacoustic tomography. Recently, it 
has been proposed to use integrating linear detectors instead of point-like transducers [6, 7, 28, 
29] . One such detector consists of a long straight segment of optical fiber that serves as a sensor 
of an optical interferometer; the measurements are proportional to the integral the acoustic 
pressure over the length of the fiber. According to the experimentalists, such detectors have 
much better sensitivity and spatial resolution than the conventional transducer. It has been 
shown [6] that the reconstruction problem in this case reduces to the inversion of a set of the 
circular mean Radon transforms followed by the inversion of a set of the regular 2-D Radon 
transforms (where the circular mean Radon transform is a set of normalized integrals over 
circles with centers lying on a curve) The same practical reasons as in 3-D case lead to the 
requirement that the centers of the integration circles lie on an open bounded curve satisfying 
(for a given ROI) the "visibility" condition. 

In the present paper we study mostly the 2-D case; in order to simplify the analysis we will 
concentrate on a truncated circular geometry as described below. The 2-D region of interest 
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centers of integration circles lie on the circular arch 'y{R^, Zright) 
R^ and zi < Zright}, where R^ > R, see Figure 1. Our goal is to reconstruct a Cq function 
f{x) supported in Q from the known values of its integrals g{z,r) over circles S(r, 2) of radii 
r centered at points z E ■y 
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(The circular mean Radon transform is defined by the normalized integrals (or means), i.e. by 
values of g{z, r)/27rr. We, however, prefer to work with the integrals g{z, r) ). The invertibility 
of the circular mean Radon transform is well known for such geometry [3]. The stability of 
the reconstruction problem is, again, determined by the "visibility" condition (see [20,24,25, 
38,39]), which for this geometry is equivalent to the inequality Zright > x right- Our goal is to 
develop an efficient computational algorithm for the solution of this problem. 



2 Outline of the method 



The present algorithm is based on precomputing the approximations of plane waves in Q by 
the single layer potentials in the form J Z{X\z — x\)p{z)dl{z), where p{z) is the density of 
the potential and Z{t) is either the Bessel function Jo{t) or the Neumann function Vo(^)- 
Specifically, given a wavevector ^ we find numerically the densities p^j{z) and p^^y{z) of the 
potentials 



Wj{x,p^,j) 



Jo{\\z ~ x\)p^^j{z)dl{z), 



WV(x,pg,y)= / Yo{\\z - x\)p^^Y{z)dl{z), 

where A = |.^|, such that 

Wj{x, p^^j) + Wy{x, p^^y) ^ exp(-i^ ■ x), \/x eQ, 



(1) 
(2) 

(3) 



Obtaining such approximations is not trivial; we discuss this issue in more detail in the next 
section. However, if the densities p^j and p^^Y have been found for all ^ then function f{x) can 
be easily reconstructed. Indeed, let us introduce convolutions Gj{\,y), Gy{^,z) as follows 



Gj{\,z) 



f{x)Jo{\\z — x\)dx, 



Gy{X,z)= / f{x)Yo{X\z-x\)dx. 
Jn 

We notice that the boundary values of these functions for all z G 7 can be computed from 
projections g{z,r): 



Gj{\,z) 



g{z,r)Jo{\r)dr, 



R^ 



Gy{\,z)= / giz,r)Yoi\r)dr. 
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(4) 
(5) 



Consider now the Fourier transform f{^) of f{x) 
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f{x) exp{—iC, ■ x)dx. 



Using (3) we obtain 
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/(x) [Wj{x, p^j) + WYix, p^,y)] dx 
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f{x)Jo{X\z — x\)dx 



f{x)Yo{\\z~x\)dx 



P(.j{z)dl{z) 
p^y{z)dl{z) 



[p^,j{z)Gj{\ z) + p5,y(2;)G'y(A, z)] dl{z). 



(6) 



Formulas (4) and (5) in combination with (6) allow us to reconstruct (approximately) the 
Fourier transform f{^) of f{x). Now f{x) can be recovered by inverting the 2D Fourier trans- 
form. In order to minimize the operation count we will compute Gj{X, z) and Gy (A, z) (using 
(4) and (5)) for a set of fixed values of A = |.^|; thus, the values of /(^) will be found on a 
polar grid in ^. One can now use a high-order 2-D interpolation to obtain values of the Fourier 
transform on the Cartesian grid and apply the inverse 2D Fast Fourier Transform (FFT). 
Another approach is to utilize the famous slice-projection theorem (see, for example [21]) and 
reconstruct from the values of /(^) on the polar grid the regular Radon projections of f{x)\ 
the function then can be obtained by the application of the FBP algorithm [14, 21]. We have 
chosen the latter approach. The whole algorithm can be briefly outlined as follows. 

1. Choose a uniform polar grid in ^. For each value of ^ij = Xi{cos6j,sin6j) precompute 
p^. .^j and p^- .y as described in the next section. This step does not depend on values 
of g{z,r) and needs to be performed only once for each particular geometry. 

2. Given values of g{z,r), compute Gj{\i,z) and Gy{K,z) (using (4) and (5)) for each 
value of i. This can be done using the trapezoid rule in r. 

3. Compute f{^ij) according to (6), using the same discretization of the integral as was 
used to obtain (1) and (2). 

4. Apply 1-D FFT to values of f{C,i,j) for each fixed j, to reconstruct the classical Radon 
projections of f{x) corresponding to the angular parameter 6j. (Alternatively, in the 
presence of strong noise in the data one can introduce a lower-pass filter 77 (^) and use 
vilCijl) f i^hj) instead of /(^j) on this step). 

5. Use the well-know FBP algorithm [14, 21] to reconstruct f{x) from the standard Radon 
projections. 

If we assume, for simplicity, that the size of the reconstruction grid is n x n, and that 
the number of detectors and the number of integrals in one projection are of the same order 
(say, 0{n)), then steps 2, 3, and 5 of the algorithm require O(n^) floating point operations; 
step 5 is faster ((9(n^logn) flops). Step 1, described in the next section, is more expensive 
computationally. The precise operation count depends on the algorithm used to compute n 
Singular Value Decompositions (SVD) for matrices of size (nxn), and can be 0{n'^) operations 
or higher. However, step 1 needs to be performed only once for each particular geometry. The 
number of data generated on step 1 and stored on a hard drive (values of the densities p^i^j 
and p^^ ^y) is of order 0{n^). The rest of the algorithm, including reading the precomputed 
densities, is completed in 0{n^) operations (similarly to the FBP). 

Obviously, the feasibility of this method hinges on our ability to find approximations in 
the form (3). In the following section we discuss the existence, accuracy and stability of such 
approximat ions . 



3 Approximations of plane waves by the single layer 
potentials 

3.1 Full circle acquisition 

Before studying the more interesting case of an open acquisition curve 7 let us consider a 
simpler case of the full-circle acquisition geometry. Suppose the 2-D region of interest Q is the 
disk of radius R centered at the origin and the detectors lie on the concentric circle 7 of radius 
R^ > R (this formally corresponds to a particular case when Xright > R and Zright > -^7)- In 
order to represent the plane wave exp{iC, ■ x) by the single layer potentials supported on 7, we 
expand the wave in the Fourier series in polar angle 6^ as follows: 



+00 



exp(i^ -x) = ^ rj|„|(A|a;|) exp(m [9^ - 9^])) 



(7) 



where C, = A(cos^^, sin^^), x = \x\{cos9x,sin9x) (this is a well known Jacobi- Anger expansion 
[9,35]). We also utilize the addition theorem [9,35] for the Hankel function Hq '{■) = Jq{-) + 

^Yoi■y. 

+00 

h!,'\\\x - z\) = Y. iJ«(A|z|)JH(A|x|)exp(m[^.-^J), (8) 

n=— 00 

where z = |2;|(cos^2, sm9z). By integrating equation (8) with the density exp(m^2) supported 
on 7 one obtains 



i^Q {\\z — x\) exp{in9z)dl{z) = 27ri?^iJ|L (Ai?^) J|„|(A|a;|) exp(m^2 



or 



'H 



H^HXRj) / H^o\M^ - ^1) exp{tn9,)dl{z) = 27tR-, H^{\R^) J|„|(A|x|) exp(m^,). 
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The latter expression can be used to express the cylindrical wave J\n\{)\\x\) exp{in9x) in the 
form 



J\n\{X\x\)exp{in9^ 




Jo{\\z — x\) exp{in9z)dl{z) 



yo(A|-2 — x\) exp{in9z)dl{z). 



Further, by substituting this formula into Jacobi- Anger expansion (7), one can represent 
exp(i^ ■ x) as the sum of the single layer potentials 



exp(i^ ■ x) 



1 



2nR 
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Wj{x,Pi:j)+WY{x,Pi:y^ 



i''exp{in[e,-e^])) 



Jo{\\z — x\)dl{z) 



Yoi\\z - x\)dl{z) 



(9) 



with the densities p^,j{z) and p^x{z) defined by the formulas 



Pui^) 



Pi,Y[z) 



2nR. 
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-R ^-^ 



J\n\{\R-y j .^ 

Yj„|(A-R. 



2" exp(2n[^, -%])), 
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The above series converge uniformly due to the fast growth of the Hankel functions (see, for 
example [9]) as the order n goes to infinity. Interestingly, in this simple case of the full-circle 
acquisition, representation (9) is exact, and therefore the reconstruction algorithm described 
in the previous section is also theoretically exact (in this particular case). 

For future use we compute the L^ norm A^(A) of the pair of the densities {p^j, P(,y) defined 
as follows 



[plA^)+plYi^)]dl{z). 
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N\X)^\\ip^j,p^,y)\\l. 

Due to the orthogonality of the complex exponents one obtains the following simple formula 

+ 00 



N\X) = J2 



1 



\^\n\ i^R-y) 



2- 



(10) 



The values of N{X) can be easily computed numerically; it turns out that for large values of 
A this function grows asymptotically linearly and A^(A) ~ v2A. 



3.2 Open curve 7: general considerations 

The case when detectors are placed on an open curve 7 is more complicated. In particular, in 
this case equation (3) cannot be made into exact equality. Indeed, the single layer potentials 
defined by equations (1), (2) are solutions of the Helmholtz equations in M^\7; these functions 
decay at infinity as (9(|a;|'-'^/^). The plane wave with |,^| = A also solves the same Helmholtz 
equation, but it does not decay. If these two solutions were made to coincide within the open 
set Q, they would also have the same behavior at infinity, which is clearly impossible. 



Not everything is lost, however. Let us consider the problem of approximating s(x) G 
L'^{dQ) by a single layer potential Wfj(i){x, p) in the form 

Whw{x,p)= f H^^\\\z-x\)p{z)dl{z), (11) 



where Hankel function H^ {\\z — x\) coincides (up to a constant factor) with the free-space 
Green's function of the Helmholtz equation 

Am + A\ = 0, 

subject to the radiation conditions at infinity. Assume additionally, that A is not an eigenvalue 
of the Laplacian on Q, with zero boundary conditions. If s{x) cannot be approximated by 
potentials (11) in the L^ sense then there exists a non-zero function t{x) G L'^{dVL) orthogonal 
to all such potentials: 



t{x) 
an 



H^^\\\z-x\)p{z)dl{z) 



dl{x) = 0, yp{z) G L^(7). 



By interchanging the order of integration we obtain 



rWi 



t{x)H^Q>{X\z-x\)dl{x) 
Uan 



di{z) = 0, M^) e L\^), 



which, in turn, implies 

T{z)= f t{x)Hl^^\\\z-x\)dl{x) = \/ze-f. 

JdQ 

Function T[z) defined by the above expression, is also a single layer potential with density 
t{x) supported on dfl. This function is real-analytic in ]R^\r2. Since it vanishes on 7 it must 
vanish on the whole circle S(r, 0) = [x \ xf + x^ = R"^^ .Thus, T{z) is the unique solution 
of the Dirichlet problem for the Helmholtz equation [9] in the exterior of the disk B(r, 0) = 
\^x\xl + X2< R^} circle satisfying zero boundary conditions on S(r, 0) and the radiation 
condition at infinity. Therefore, T{z) identically vanishes in ]R^\]B(r, 0). By the analyticity of 
the solutions of the Helmholtz equation, T{z) vanishes in M^\f2. On the other hand, by the 
continuity of the single layer potentials, T{z) also approaches zero when z and approaches 
dQ from inside of Q. Since A is not an eigenvalue of the Laplacian on Q, T{z) vanishes in Q. 
Now the well-known jump condition for single layer potentials [9] implies that t{x) identically 
equals zero on dQ. This contradiction proves the following 

Theorem 1 An arbitrary s{x) G L'^{dQ) can be approximated by potentials in the form (11) 
in the L^ sense. (In other words, single layer potentials (11) are dense in L'^{dVL).} 

A very similar statement can be proven if one replaces Hankel's function H^ hy Hq in 
the above proof. Therefore, the sums of potentials 

WHw{x,pi) + WHi2,{x,p2)= j Hi^\\\z-x\)pi{z)dl{z)+ I Hf\\\z-x\)p2{z)dl{z) 



with arbitrary densities pi{z) and P2iz) are also dense in L'^(dfl). Alternatively, due to the 
equations 

H!,'''\t) = Mt)±tYo{t), 

the sum Wj{x,pj) + Wy{x,py) can be used for approximation, where 



Wj{x,pj)= / Jo{X\z-x\)pj{z)dl{z), 
Wy{x,py)= [Yo{X\z-x\)pY{z)dl{z), 



7 



again, under the assumption that A is not an eigenvalue of the Dirichlet Laplacian. 

Our goal is to approximate the plane waves in Q by the single layer potentials. Due to the 
uniqueness and stability of the Dirichlet problem for the Helmholtz equation, it is enough to 
approximate the boundary values of these functions on dQ, which can be done as explained 
above — except for the case when A is the eigenvalue of the Dirichlet Laplacian. In the 
latter case, one can approximate the normal derivative of the target function by the normal 
derivative of the single layer potential. A derivation similar to the one in the beginning of this 
section shows that the normal derivatives of the layer potentials ^4^ {Wj{x, pj) + Wy{x, Py)) 
also form a dense set in L^((9r2), if A is not an eigenvalue of the Neumann Laplacian on Vt. 

Finally, since the eigenvalues of Laplacian on VL may not be known in advance, it is ad- 
vantageous to use simultaneously both Neumann and Dirichlet data. Namely, in order to 
approximate a solution u{x) of the Helmholtz equation on Vt by the single layer potentials, we 
form a vector function ('u(x), Qnix) '^^^))\an ^"^^ ^^J" ^^ approximate it in the L?' sense by the 

functions in the form \Wj{x,pj) + Wy{x,py), q^ \Wj{x,pj) + Wy{x,py)\ 



an 



3.3 Stability and regularization 

As established in the previous section, given a plane wave exp(— i^ ■ x), one can can find a 
sequence of densities [p\ j, p\y), k = 1,2,..., such that the sum of potentials W,j{x, p^ j) + 

Wy{x^ pi y) converges to the wave as described. However, in the case of open 7 the sequence of 
densities themselves cannot have a limit in L^ sense; if this limit existed, the sum of potentials 
would exactly equal the plane wave, which cannot happen, as discussed in the beginning of 
section 3.2. 

Moreover, the sequence [P^ j^ PfYJ can start growing very fast in the L^ norm defined by 
the formula 

\\iPJ,PY)\\l= [[\pj{z)f + \pY{z)f]dl{z) 

This, in particular, should occur in the situation when curve 7 does not satisfy the visibility 
condition. Indeed, it is known that the reconstruction problem is strongly unstable, if this 
condition is not satisfied [20,24,25,38,39]. On the other hand, possible instabilities in the 
present algorithm are associated with densities that are strongly oscillating and large in the L^ 
norm. Observe that in equation (6) we compute the inner product of densities p^j, p^^y with 
functions G'j(A, z) obtained from the measured data. Any component of noise well correlated 
with p^^j and p^y will be strongly amplified if these functions are large (in L^ norm). This 
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amplification is the only source of instability in reconstructing the value of f($,) at a particular 
point ^ by our method, and it has to occur to render the reconstruction unstable — as it should 
be in accordance with the theory. 

Our goal, of course, is to use the present technique with the data acquisition configurations 
that satisfy the visibility conditions. It is difficult to find theoretically a sharp estimate on the 
behavior of some "optimal" sequence of densities that would combine convergence of potentials 
with the slow or no growth of the densities. Instead, we propose a regularized algorithm for 
computation of "good" approximations. 

Let us introduce the families F and T of pairs of L^ functions defined on 7 and dQ: 

^ = {{qiiz),q2iz))\qi,q2 e L\^)}, 

r = {{p,{z),p2{z))\p„p,eL\dQ)}. 

Define the inner products for functions from F and T as follows 



Vq,sGF, 



(q,s)r= / 



qi[z)si[z) + q2[z)S2[z) 



dl{z), 



Vp,reT, (p,r)^ 



an- 



Pi{x)ri{x) + p2{x)r2{x) dl{x). 



Further, let us introduce the operator A that maps a pair q gF into a pair p gT according 
to the following formula 



p = A(q) = [Wj{x, qi) + Wy{x, gs), 



1 d 



A dn{x 



[Wj{x,qi) + WY{x,q2)] 



dQ 



where single layer potentials Wj{x, gi), Wy{x, ^2) are defined, as before, by equations (1), (2). 

Given boundary values of the plane wave u^{x) = (exp(i^ ■ 
like to find a pair of functions (densities) {p^^j{z), p^^y{z)) such that 



^^ exp(i^ ■x))an we would 



A((p5,j,P5,y)) ^u^. 



(12) 



subject to the requirement that the densities are "not too large". There is more than one 
way to find such regularized solutions; we utilize the SVD of the operator A. Namely, we find 
(numerically) the sets of pairs (left and right singular vectors) q'--'-* G F, p*^-'^ G T, j = 1,2, .., 
and singular values aj, ai > (72 > ... > a-n > ..., such that 

(q»,q(^-)>^ = <5,„ 
(p»,pO-))^ = 5^^^., 

P^'^ = ^^■A(qj), 

where 5ij is the Kronecker symbol {5i^j = if i 7^ j, and 5i^j = 1 otherwise). Now the desired 
approximation is given by the sum 



Jmax 



(Pc,jW,Pc,yW) ~ X1^^^^(^)~(^«'P^^^)t' 



(13) 



j=i 
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where jmax serves as a regularization parameter. The increase in j'^ax yields a closer L^ fit of 
the single layer potential to the plane wave, but it also may (and in certain cases will) lead 
to the unbounded growth of the densities. 

A frequently used regularization technique for solving ill-posed problems using the SVD 
decomposition is to drop from (13) all the terms with ctj smaller than a certain threshold 
(Tmin (and thus define the jmax)- However, it is not clear how to choose such (Tmin optimally. 
Moreover, in general (Tmin should depend on the frequency A, since the norm of the densities 
may grow with A as we saw in the example of circular acquisition geometry. 

Instead of chosing (Tmin, we propose to use the circular case as a benchmark to determine the 
regularization parameter jmax- We found in the numerical experiments that very good results 
are obtained when jmax is chosen to be the largest number such that 1 1 (p^,j, p^y) 1 12 < KN{\), 
where A^(A) is computed using equation (10) for the circular case with the same values of R 
and R-y, and i^ > 1 is a constant; in all numerical experiments presented below K was equal 
to 1.5. 

The most computationally expensive step in finding the densities is the SVD decomposi- 
tion; however, the same SVD is used for all plane waves with the same frequency A = |.^|. 
Algorithmically, we computed the SVD by first discretizing the integrals that define operator 
A, and by enforcing equation (12) at some set of collocation points on dVt. This results in a 
matrix that represents a discretized version of A; the SVD is then computed using subroutine 
DGESVD from LAPACK. In practice, the discretization of the integrals is dictated by the 
number and location of the detectors; we assumed that they were distributed uniformly over 
7. We also utilized equispaced collocation points on dVL\ the number of points was chosen to 
be twice the number of the detectors. 

One can notice that although this approach guarantees bounded solutions {p^j,p^y) , 
there is no theoretical estimate on the accuracy of the approximation of the plane waves by 
the corresponding potentials. However, when the densities have been found, it is very easy 
to compute numerically the approximation error, and, if necessary, to re-run the computation 
with different values of parameters. 

The computation of the densities {p^j{z), p^y{z)) constitutes the first step of the recon- 
struction algorithm (see section 2). It is rather time consuming; however, for a particular 
acquisition geometry it needs to be done only once. 

4 Numerical examples 

In order to verify the accuracy and stability of the present algorithm we conducted a series 
of numerical experiments . In particular, we studied two acquisition geometries (which we 
will call #1 and #2) corresponding to different values of parameters Xnght and Zright, (see 
Figure 1). In both geometries R and R^ are equal to 1 and 1.3 respectively. Geometry ^1 
has Xright=Zright = 1; in other words, the ROI is a unit circle. Geometry #2 is defined by 
values Xnght=Zright = 0, which corresponds to a half-circle acquisition curve and a half-disk 
ROI. Both geometries satisfy the visibility condition. In all experiments reconstruction was 
performed on an equispaced Cartesian 129 x 129 grid, from simulated projections corresponding 
to 500 equispaced detectors, each measuring 129 circular integrals with equispaced radii. (The 
radial step in the projections coincides with the step of the reconstruction grid under such 
discretization). 
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Figure 2: Geometry ^1: (a) phantom, (b) reconstruction from the accurate projections (c) 
effect of 15% noise in projections (d) reconstruction from the noisy projections with the 
additional filter 

The goal of our first test was to verify the feasibility of the approximation of plane waves 
by the single layer potentials. We found that the most difficult plane waves to approximate 
were the ones propagating in the vertical direction. In geometry ^1 the method described 
in Section 3.3 with K = 1.5 gave approximations accurate up to 4 decimal places or better. 
In particular, for the plane wave propagating in the vertical direction with the wavelength 
corresponding to the Nyquist frequency of our 129 x 129 grid, the maximum pointwise error 
did not exceed 8 ■ 10^^. Similar (although slightly less accurate) results were obtained in the 
geometry ^2. 

These accurate approximations of the plane waves lead to very accurate reconstructions 
of smooth images. For geometry ^1 we defined a smooth phantom as follows. First, we 
introduced a Cq(]R) function h{t): 



h{t) 



Co lo 




-1*1 



sm 



(7rs)(is 



< |t| < 1 

Itl > 1 



where constant Cq was chosen so that h{0) = 1. Function h{t) is even, compactly supported on 
[—1, 1], 8 times continuously differentiable on M function, with h{l/2) = 1/2. We then used 
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this function to construct the 8 times continuously differentiable phantom f(x) = h ( ^ ^^ ) + 

h ( ~^ 1 , as a sum of two bell-shaped rotationally invariant functions. For geometry ^1 

we used parameters Xi = (0.3, 0.3),a;2 = (—0.4,0.2), ri = 0.55, r2 = 0.5. Reconstruction 
using the present method (without additional filtration) resulted in maximal point-wise error 
of 7.3 ■ 10^^. Similarly small reconstruction errors were obtained for smooth phantoms in 
geometry #2. We do not present here the gray scale pictures of these phantoms and the 
corresponding reconstructed images since they look identically. The experiments with smooth 
phantoms show that our method indeed reconstructs accurately all the spatial frequencies of 
the image (including the lower ones) — as opposed to parametrix approximations. 

Our next experiment was with discontinuous phantoms. Point-wise accurate reconstruc- 
tions are not possible with such functions due to the aliasing errors and the Gibbs phenomenon; 
the goal is to obtain images that appear qualitatively correct, with low noise amplification. As 
a phantom in geometry ^1, we used the sum of the characteristic functions of circles as shown 
in Figure 2(a). Part (b) of the latter figure demonstrates the image reconstructed from accu- 
rate projections. In order to analyze the sensitivity of the method to non-exact measurements 
we added to the projections white noise with intensity of 15% of the signal (in L^ norm); the 
reconstruction is shown in Figure 2(c). As it is frequently done in tomography, in order to 
reduce the effects of noise one can apply a low-pass filter on step 4 of the algorithm. Part 
(d) demonstrates the effect of such additional filtration, with filter ?7(^) = cos{^\^\ / \Myquist) , 
where \Nyquist is the Nyquist frequency of the reconstruction grid. 

Quite similar results were obtained in geometry ^2. As a phantom we used again a sum 
of the characteristic functions of circles supported within the ROI (the left half of the unit 
disk)), as shown in Figure 3(a). Parts (b), (c), and (d) of that figure demonstrate, respectively, 
reconstruction results from the accurate projections, noisy projections, and reconstruction 
from noisy projections with additional filtration (using the same filter as in the previous 
example). 

In order to compare stability of our algorithm to that of the classical FBP, we have recon- 
structed the two phantoms (shown in Figure 2(a) and Figure 3(a)) using the latter method, 
from a similar number of the standard Radon projections with the same level of noise (not 
shown here). The general quality and the level of noise in the reconstructions was quite close 
to those of the images obtained by the present method. 

From the practical point of view it is interesting to know how the reconstruction is affected 
by acoustic sources located outside of the ROI. It is known that, with the exception of the 
methods based on expansion in the eigenvalues of the Dirichlet Laplacian on a closed domain [2, 
19], all the other exact reconstruction techniques will produce incorrect results in the presence 
of such sources (see [1, 13, 16] for further discussion of this phenomenon). In our case, such 
sources will be present if the electromagnetic wave impinges on the parts of the patient's 
body located outside the ROI. In other to model this situation in geometry ^^2 we used the 
previously used phantom supported within the unit disk, as shown in Figure 2(a). The result 
of the reconstruction in the absence of noise is shown in Figure 4(a); part (b) demonstrates 
the effect of 15% noise and the additional filtration. One can notice in this figure that the 
parts of the source located outside ROI are not reconstructed correctly (not unexpectedly). 
Nevertheless, by covering the right half of the image, it is easy to see that the reconstruction 
within the ROI is very little (if at all) affected by the presence of additional sources outside 
the region. Interestingly, in the right half of the image (outside ROI), the "visible" material 
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Figure 3: Geometry ^2: (a) phantom, (b) reconstruction from the accurate projections (c) 
effect of 15% noise in projections (d) reconstruction from the noisy projections with the 
additional filter 

interfaces are reconstructed quite well, while the "invisible" ones (such that the normal does 
not intersect the measuring surface) are noticeably smeared. 

In accordance with the theoretical operation count, step ^1 turned out to be the most 
expensive part of the algorithm; in the experiments described in this section the computation 
of the densities took several hours. However, this step needs to be performed only once for a 
given geometry. On consecutive runs of the reconstruction program, once the pre-computed 
densities had been read from the hard drive, the algorithm took a fraction of a second to 
complete. The time required to read the densities was about 4 seconds. If our method were 
to be used to process the data obtained by an acquisition scheme with linear detectors (which 
requires the inversion of a set of 2-D problems), the densities would need to be read only once, 
and each of the 2-D problems would be inverted in a fraction of second. 

5 Concluding remarks 

We have presented an efficient reconstruction algorithm applicable to problems of thermoa- 
coustic tomography with the detectors lying on an open curve. The method is based on L^ 
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Figure 4: Geometry ^2: (a) reconstruction with some sources outside of the ROI (compare 
to Figure 2(a)) (b) the same with 15% noise in the projections and with the additional fiher 

approximations of plane waves by certain single layer potentials; we have proven that such 
approximations are possible if the measurement curve is an open circular arch. We have also 
verified numerically that, for the truncated circular geometry satisfying the "visibility" con- 
dition one can obtain very accurate approximations with bounded densities, which, in turn, 
leads to stable image reconstructions. 

In conclusion, we would like to add several remarks: 

• The theorem we have presented does not restrict the shape of the ROI; it can be arbitrary. 
Of course, for stable reconstruction the combination of the ROI and the acquisition 
surface should satisfy the "visibility" condition. 

• In the proof of the theorem we used the fact that if the solution of the Helmholtz 
equation vanishes on a circular arch, it must vanish on the whole circle, and, therefore, 
in the exterior of the disk. A similar statement is also true if the acquisition curve is a 
part of any closed analytic curve, and the theorem extends to such configurations. The 
theorem also remains valid if the measurement curve is a continuous non-analytic curve. 

• In order to use this approach with non-circular acquisition curves, the regularization 
technique we presented may require some modifications, since it is based on the com- 
parison with the full-circle case. 

• The technique we propose can also be used in 3-D, for image reconstruction from the 
detectors lying on an open surface. Compliance with the "visibility" condition is, again, 
necessary for a stable reconstruction in this case. We anticipate that the pre-computation 
of the densities may become quite time-consuming in 3-D, due to the increased dimen- 
sionality of the problem. The development of a fast algorithm for this part of our method 
is required to make this approach practical (in 3-D). This will be the object of our future 
research. 
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